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Abstract: We developed an efficient numerical algorithm for computing the spectrum 
of anomalous dimensions of the planar A/" = 4 Super-Yang-Mills at finite coupling. The 
method is based on the Quantum Spectral Curve formalism. In contrast to Thermodynamic 
Bethe Ansatz, worked out only for some very special operators, this method is applicable 
for generic states/operators and is much faster and more precise due to its Q-quadratic 
convergence rate. 

To demonstrate the method we evaluate the dimensions A of twist operators in s[(2) 
sector directly for any value of the spin S including non-integer values. In particular, we 
compute the BFKL pomeron intercept in a wide range of the ’t Hooft coupling constant 
with up to 20 significant figures precision, confirming two previously known from the per¬ 
turbation theory orders and giving prediction for several new coefficients. Furthermore, we 
explore numerically a rich branch cut structure for complexified spin S. 
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1. Introduction 

Many years of exploring integrable structures in planar AA = 4 super Yang-Mills and its 
AdSs X S® string dual have led to a remarkably simple system of equations for the exact 
spectrum of the theory known as Quantum Spectral Curve (QSC) equations. They are 
expected to capture the conformal dimensions/string state energies at any value of the ’t 
Hooft coupling [1, 2]. The QSC equations are formulated as a set of Riemann-Hilbert type 
equations for a few functions. As a limiting case this system incorporates the renowned 
asymptotic Bethe ansatz, but also includes all wrapping corrections essential for finite 
length operators and perhaps constitutes the ultimate solution of the spectral problem. 
The Quantum Spectral Curve has a transparent algebraic origin related to the underlying 
psu(2, 2|4) symmetry of the problem. It was derived in [1, 2] from the discrete integrability 
of Hirota dynamics underlying the Y-system enhanced with a specific analyticity condition 
worked out in [3, 4, 5, 6, 7]. In contrast to some explicitly known integral forms of the 
Y-system [5, 8, 9, 10, 11] the QSC is applicable for any local operator or string state. 

The striking simplicity of the QSC formulation has already led to a rapidly growing 
body of exact as well as perturbative results. First, in [1, 12] it was pointed out that it can 
be solved analytically in a near-BPS regime, e.g. for expansion in the coupling A or in the 
spin S, in particular in [12] the anomalous dimensions for twist operators in the 51(2) sector 
were computed at any coupling to order S'^, giving new analytical predictions at strong 
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coupling for some local operators and the BFKL pomeron intercept. As an extension of 
the observation of [1] powerful techniques for expansion in the coupling were developed in 
[13, 14] which led, impressively, to 10-loop anomalous dimensions for numerous operators, 
as well as 6-loop predictions at any S for twist-2 operators. Finally, in [15] the analytic 
continuation to the BFKL region iS ~ — 1 was explored, and the leading order BFKL 
equation in SYM was derived in this approach. It is also important to mention that the 
QSC method was also applied to ABJM theory [16] opening a way for various explicit 
calculations in particular to that of the interpolation function h{X) [17], the mysterious 
extra component of the spectral problem in this model. 

Even though many explicit analytic results are now available both at strong and weak 
coupling, one important range of applications of the QSC that has remained unexplored 
till now is the numerical investigation of the spectrum at finite coupling. 

Previous numerical methods based on TBA, even limited to a few operators^, low 
precision and slow convergence rate gave, nevertheless, several highly important results, 
allowing, in particular, the computation of the anomalous dimension of a nonprotected 
(Konishi) operator in a planar 4d theory at finite coupling [18]. Numerics also gave a 
prediction for the strong coupling Konishi anomalous dimension which was later confirmed 
by several methods [19, 20, 21, 22, 12, 23, 24, 25] . The main goal of the present work is 
to remove the limitations of the previously known methods by developing an algorithm for 
a numerical solution of the QSC. 

The low precision and performance of the TBA-like approach was mainly due to the 
complicated infinite system of eqnations and cumbersome integration kernels. The QSC 
includes only a few unknown functions and thus can be expected to give highly precise 
numerical results. However, the QSC equations are functional equations supplemented 
with some analyticity constraints of a novel type which makes it a priori not a trivial task 
to develop a robust numerical approach. 

In this paper we propose an efficient method to solve the QSC numerically and illustrate 
our method by a few examples. Among the several equivalent formulations of the QSC we 
identified the equations which are best-suited for numerical solution^. We implemented onr 
algorithm in Mathematica and were able to get a massive increase in efficiency compared 
to the TBA or FiNLIE systems [18, 8, 24, 9]. With one iteration taking about 2 seconds 
we only need 2 — 3 iterations (depending on the starting points) to reach at least 10 digits 
of precision. Quite expectedly, the precision gets lost for very large values of the’t Hooft 
coupling. Nevertheless, without any extra effort we reached A ~ 1000 keeping a good 
precision, which should be more than enough for any practical goal. 

Not only does our approach work for any finite length single trace operator and in 
particular for any value of the spin, it also works with minimal changes even away from 
integer qnantum nnmbers! We demonstrate this in the particularly interesting case of 
the s[(2) twist-2 operators. Their anomalous dimension analytically continued to complex 

^only for a few operators the complicated structure of the “driving terms” was deduced explicitly in a 
closed form. Even for those operators the driving terms may change depending on the value of the coupling. 

^one may call this sub-system of equations as Pcj-system, in contrast to previously used P/r-system or 
Qcu-system 
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Figure 1: Riemann surface of the function S{A) for twist-2 operators. Plot of the real 
part of S'(A) for complex values of A, generated from about 2200 numerical data points for A « 6.3. 
We have mapped two Riemann sheets of this function. The thick red lines show the position of 
cuts. The upper sheet corresponds to physical values of the spin. Going through a cut we arrive at 
another sheet containing yet more cuts. 


values of the spin S is known to have a very rich structure, in particular the region S' ~ — 1 
is described by BFKL physics. As we show, within the framework of QSC it is not hard 
to specify any value of the Lorentz spin S as the conserved charges enter the equations 
through the asymptotics which can in principle take any complex values. Then we can 
compute the analytically continued scaling dimension A directly for complex S (or even 
interchange their roles and study S as a function of A). The result of this calculation can 
be seen on Fig. 1. 

Let us stress that the algorithm is very simple and mainly consists of elementary matrix 
operations. As such it can be easily implemented on various platforms. In particular, we 
believe the performance could be increased by a few orders with a lower level, e.g. C++, 
implementation. In this paper we mostly aim to demonstrate our algorithm, prototyped 
in Mathematica. To illustrate how the algorithm works, as an attachement to this paper 
we provide a Mathematica notebook with a simple implementation of our method. 

Finally, to improve the performance of our method we used the further simplification 
of the QSC obtained in [26], which allows us to eliminate auxiliary functions coij from our 
algorithm and close the equations using Q-functions only (we demonstrate this for the sl{2) 
sector states). 
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The structure of the paper is as follows. In section 2 we introduce the Quantum 
Spectral Curve and discuss the key equations for the numerical implementation. In section 
3 we describe how our algorithm works in general. We also demonstrate our method on 
specific examples for twist-2 operators in the s[(2) sector, starting with physical operators 
(such as the well studied Konishi operator). Then in section 4 we explain how to extend 
the method to non-integer values of the spin. We discuss Fig. 1 in a bit more detail, and 
also describe our high precision evaluation of the BFKL pomeron intercept^. The final 
section contains our conclusions and speculations. 

2. Review of the AdS/CFT Quantum Spectral Curve 

Let us start by introducing the main equations of the Quantum Spectral Curve framework. 
For a full review of the QSC we refer the reader to [2], but we will try to make the paper 
self-contained. 

For integrable spin chains Baxter Q-polynomials play a central role. They carry com¬ 
plete information about the state and the spectrum and can also be thought of as wave 
functions of the model. They are completely determined by a set of functional relations 
(known as Q-system) and by the polynomiality condition. 

In AA = 4 SYM the situation is similar. Its spectrum of anomalous dimensions is 
described by a set of Q-functions satisfying the same Q-system relations. However, in this 
model the Q-functions are no longer polynomials, but are analytic functions of the spectral 
parameter u with branch cuts. The positions of the branch points depend on the ‘t Hooft 
coupling A: they are situated at ±2^ -|- in, n £ Z, where g = ^- In the limit when A is 
small the system reduces to the usual psu(2,2|4) Heisenberg spin chain. For the class of 
functions with cuts the Q-system alone is not constraining enough and one has to specify 
the monodromies of the Q-functions in order to close the system. These monodromies 
were found in [I] and take the form of a Riemann-Hilbert problem as described below. The 
resulting set of equations found in [1, 2] is known under the name Quantum Spectral Curve 
(QSC). 

The Q-system of A/" = 4 SYM is composed of 2® Q-functions. However, the algebraic 
relations between them allow us to choose a much smaller subset, which will be complete in 
the sense that the rest of Q-functions can be generated from the selected ones algebraically. 
A convenient choice for such a subset consists of 4 -|- 4 functions Pa (it) and Cli{u) {a,i = 
1,..., 4). One can say that Pa describe the degrees of freedom whereas Qj correspond 
to the AdS^ part. A particularly nice property of P’s is that they have only two branch 
points at ±2g when they are connected by a “short” cut [—2g] 2g\ (see Fig. 2). This means 
that there are more branch points on the next sheet, but for this choice of the cut they do 
not appear on the first sheet. Very similarly Q’s have only two branch point on the main 
sheet if the cut is taken to go through infinity. In a sense this reflects the non-compactness 
of the AdS^ part of the space. 

Whereas the coupling determines the position of the branch points, the quantum num¬ 
bers of the state are specified through the large u asymptotics of Q-functions. Pa encode the 

®i.e. the value of S for which S{A) — 0 
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Figure 2: and Qi have one cut on the real axis in the representations with short and long cuts 

respectively. The ellipse shows the region of convergence of the series (2.4). 


compact bosonic subgroup SO{G) quantum numbers (Ji, J 2 , J 3 ), while Qi give the 50(4, 2 ) 
charges (A, 5 i, 52 ), which include the conformal dimension of the state A. Explicitly 






(2.1) 
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Note that for the numerical implementation P^ are more handy, in the sense that they 


can be expressed as a series in the Zhukovsky variable x{u) defined by u = g{x + 1/®), 

00 

p„(u) = y yy . (2.4) 

^ x^(u) 
n=Ma 

This series is convergent everywhere on the upper sheet and also in an elliptic region around 
the cut on the next sheet (see Fig. 2). A similar parametrization for Qj will not cover even 
the upper sheet. Fortunately, in the whole set of 2 ® Q-functions there are other 4 functions 
with one single cut, which are denoted as P“(m), a = 1,... ,4. Together with Pa(rt) they 
also form a complete set of Q-functions. In particular, one can reconstruct Qi from them. 
The procedure for this, which will be crucially important in our numerical implementation, 
is the following: 

• Find a set of 16 functions Qa\i, satisfying 

Qa\i{u+^^)-Qa\i{u-\) = -¥a{u)¥\u)Q^,{u+h) ' ( 2 - 5 ) 

Note that this is a 4-th order finite difference equation, which entangles all Qa\i with 
fixed i. Different values of i label the 4 linearly independent solutions of this equation. 
One could also equivalently use Qh\i{u — |) in place of Qb\i{u + |) in the r.h.s., due 
to the constraint [ 2 ] 

PaP“ = 0 . (2.6) 
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This constraint also fixes some of the coefficients Ca,n- 

• The matrix can then be used to pass to Q* from P“’s, 

Q,(u) = -P“(u) Q,|i(u + i/2) . (2.7) 

The equations (2.5) and (2.7) are simply two of the Q-system relations as explained 
in [2]. We also introduce a matrix such that = —Sj and use it to define 

Q’s with an upper index: 

Q*(u) = +Pa(n) Q“l*(u + i/2) . (2.8) 

Note that since Qa\iiu) is analytic in the upper-half-plane we can also analytically 
continue these relations around the branch point at tt = 2g' to get 

Qi{u) = -P“(u) Qa\i{u + i/2) (2.9) 

Q*(u) = +Pa(u) + i/2) (2.10) 

where the tilde denotes analytic continuation to the next sheet. 

Let us also mention that large u asymptotics of P“ and Qj read [2] 

pa _ Qi _ _ (2.11) 


Since Qj can now be recovered from P^ and P“ it is not surprising that actually all 
information we need, in particular all the charges (including those in AdS^), are encoded 
in P’s alone, through 


.-Ma 


A^u 


Ma-l 


7l“°A 


ao 


^ Uj{Mao-Mj) 

'Ub^aoiMao-M,)' 


( 2 . 12 ) 


where Mj and Ma are defined in (2.2), (2.3) and uq takes values 1,2, 3,4 but there is no 
summation over oq in l.h.s. In particular, one can extract A from the last equation. 

The coefficients Ca,n and corresponding coefficients c®’” of the expansion of P“(u) need 
to be found. The constraint (2.6) fixes some of them (for example, we can use it to fix all 
ci^n)- The condition (2.12) gives the leading coefficients ^ • The remaining coefficients 
should be fixed from the analyticity constraints on P’s as prescribed by QSC. Let us 
describe these constraints. The analytic continuation of Pq to the second sheet, which we 
denote by Pq, in terms of our ansatz (2.4) becomes simply 


Pa(w) = ^ Ca,nX^{u) . (2.13) 

n=Ma 


According to [2] we should have 


Pa(^^) = i^abiu)P^ , P'^iu) = M“'(u)Pfe(u) (2.14) 
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Figure 3: ^ab is periodic as a function with long cuts, and ujij as a function with short cuts. 

where ^-ab{u) is an antisymmetric matrix with unit Pfafhan, i-periodic as a function with 
long cuts, with the discontinuity fixed in terms of Pq 

lji-ab{u) - ^iah{u) =^aPh-'PaPb ■ ( 2 - 15 ) 

Knowing the r.h.s. of (2.15) it is straightforward to reconstruct ^ab itself using the 
spectral representation of Hab 

^^ab{u) = I 

In theory one could reconstruct fiab for some fixed coefficients Ca,n from (2.15) and then 
impose (2.14) to fix the unknown Ca,n- However, in practice that is very hard to do as the 
power series (2.4) is only convergent inside the region shown on Fig 2. The problem is that 
this region does not cover the entire cut of //, which stretches to infinity. In other words, it 
would be very hard to reconstruct fiab from some given coefficients Ca,n in this direct way. 

We found that it is much more advantageous to close the analyticity conditions at 
the level of Qj, which obey very similar equations. The rule is quite simple - one has to 
interchange short and long cuts. That is, we have to introduce an i-periodic with short 
cuts function uJij (u) such that 

Q^ = coiJQJ' , Q'=a;'^Q,- , (2.17) 

^ij — QiQj QiQj • (2*18) 

with Pf w = 1. Now we see that to recover cjij one only needs to know its discontinuity 
on the interval [—2g,2g\, which is completely inside the region of convergence on Fig. 2! 
Thus the discontinuity of u can be expressed in terms of the coefficients Ca,n via (2.7) and 
(2.9), provided we know how to solve (2.5) for arbitrary and P“. In the next section 
we will describe an algorithm which allows to solve (2.5) very efficiently and then find the 
coefficients Ca,n) which yields the solution of the QSC. We will also show that actually 
finding ojij is not necessary and we can close the system in terms of just Pq, Q* and Qa\i, 
thus speeding up the calculations. 


p -2g 

dv + dv 


2q 


coth(7r(tt-'i;)) Pa('(;)Pfe('(;) - Pa(u)Pb(?;) Tperiodic, (2.16) 
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3. Description of the Method 


3.1 Step 1: Solving the equation for 

As we explained above the quantity Qa\i is at the heart of our procedure. In this section 
we will demonstrate how this set of 16 functions can be found for arbitrary and P“. 
In this procedure the precise ansatz for P is not important. However, as we will see later, 
we should be able to compute Pa(u)P^(?x) on the upper sheet for u with large imaginary 
part. Of course, having the ansatz in the form of a (truncated) series expansion (2.4) we 
can easily evaluate it everywhere on the upper sheet numerically very fast. 

The process of finding Qa\i is divided into two parts. Firstly, we find a good ap¬ 
proximation for Qa|j at some u with large imaginary part (in the examples we will need 
Im u ~ 10 — 100). At the next step we apply to this large u approximation of Qa\i a 
recursive procedure which produces Qa\i at u ~ 1. 


Large u approximation. For Im u ~ 10 — 100 we can build the solution of (2.5) as 
a Ifu expansion. This is done by plugging the (asymptotic) series expansion of Qa\i into 


(2.5): 


Qa\i{u) = U 




N 

E 

72=0 


B. 


a|2,72 


U" 


(3.1) 


where N is some cutoff (usually ~ 10 — 20). This produces a simple linear problem for 
the coefficients which can be even solved analytically to a rather high order. The 

leading order coefficients of Qa\i can be chosen arbitrarily. After that the linear system of 
equations becomes non-homogenous and gives a unique solution in a generic case.'^ 


Finite u approximation. Once we have a good approximation at large u we can simply 
use the equation (2.5) to recursively decrease u. Indeed defining a 4 x 4 matrix 


UaHu)=S'f + I>a{u)P\u) 


(3.2) 


we have 

Q„|,(u-i) = [/,'(u)Qb|i(u + |). (3.3) 

Iterating this equation we get, in matrix notation 

Qa\i{u -l) = [U{u)U{u + i)...u{u + iN)l ^ Qfc|,(u + iN+^). (3.4) 

For large enough N we can use the large u approximation (3.1) for in the r.h.s. As a 
result we obtain the functions Qa\i for finite u with high precision. 

^The matrix of this system may become non-invertible unless some constraint (which is not hard to find) 
on the coefficients Ca,n is satisfied. This constraint is fnlfilled trivially for the even in the rapidity plane 
operators considered in the next section. There is also no such problem for the situation with generic twists 
(similar to /3— or y-deformations, see the review [29]). Adding twists should correspond [2] to allowing 
exponential factors in the asymptotics of Pa and Qi, making everything less degenerate and 

providing a useful regularization. 
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3.2 Step 2: Recovering ujij 

Now when we have a good numerical approximation for Qa\i{u) we can compute Qi and 
Qj which through the discontinuity relation (2.18) will yield us LUij. 

Let us also note that, as it was argued in our paper [26], one can in fact close the QSC 
equations without calculating cuij (this was shown explicitly for the symmetric sl(2) sector 
states). This makes it possible to further speed up our numerical procedure as we will 
describe in detail in section 3.4. In the current section for completeness we will present the 
procedure without this shortcut, as for some applications it could turn out to be useful as 
well. 

One can recover uJij from its discontinuity (2.18) modulo an analytic function, as due 
to i-periodicity of uiij one can write its spectral representation as 


2g 




dvcoth(7r(u — v)) Qi(v)Qj(v) — Qi(v)Qj(v) +cu^j(u) (3.5) 


-23 


where the “zero mode” w? (u) is the analytic part of uJij — it has to be periodic, anti¬ 
symmetric in i,j and should not have cuts. We will hx it below. One can check directly 
that the term with the integral in the above expression is periodic and moreover due to 
the pole of the coth function it has precisely the correct discontinuity on the real axis to 
satisfy (2.18). Therefore this term gives a particular i-periodic solution of (2.18), while the 
io^j{u) is the general solution of the same equation (2.18) with r.h.s. set to zero, so that 
(3.5) indeed gives the general solution of (2.18). 

We note that we only need to know values of Q and Q on the cut. In our implemen¬ 
tation we use a finite number of sampling points on the cut given by zeros of Chebyshev 
polynomials. One can then fit the values of QiQj — QiQj at those points with a polynomial 
times the square root After that we can use precomputed integrals of the form 

J- 2 g cotll{^^{ui — v))v^^yv‘^ — Ag'^dv to evaluate (3.5) with high precision by a simple matrix 
multiplication, which produces the result at the sampling points ua in an instant. 

One more point to mention here is that in our implementation we only compute = 
— ojij) at the sampling points to avoid the problem of dealing with the singularity 
of the integration kernel. Note that cu’J® can be used instead of Uij in the equations like 
(2.17), because the difference is proportional to QjQ* which is zero similarly to (2.6), as 
can be shown by combining (2.6) with (2.7), (2.8). 

Finding zero modes. It only remains to fix Lo^j{u). First we notice that for all physical 
operators Uij should not grow faster than constant at infinity [2]. As the integral part 
of (3.5) does not grow either and since uj^j{u) is Aperiodic it can only be a constant. To 
fix this constant we use the following observation [2]: the constant matrix afj which Uij 
approaches at u —)• -foo and the constant matrix a~j which it reaches at u —)■ — oo are 
restricted by the quantum numbers [2] . To see this we can pick some point on the real axis 
far away from the origin and shift it slightly up into the complex plane, then from (2.17) 
we have 

ujijQ^{u + iO) = aljQ,^{u + fO) = Qi{u + iO) = Qi(u - iO). (3.6) 
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Similarly for —u we get 


(3.7) 


a^jQ^{-u + iO) = Qi{-u - iO). 

Next, notice that since Q-^ is analytic everywhere except the cut on the real axis, it can 
be replaced by its asymptotics above the real axis, i.e. Q-^(u + iO) ~ , and also 

Q'^ (—tt + iO) ~ , as we find from the previous expression by a rotation by 

TT in the complex plane. As a result we get the asymptotics of Qj at infinities and slightly 
below the real axis 


Qi(n - iO) = , Qi(-w - iO) = a-^B^u-^^. (3.8) 


Using that they are related by the analytic continuation in the lower half plane the hrst 
equation also gives 


Qi{-u - iO) = . 


(3.9) 


Combining this with (3.8) we get a relation between the constant asymptotics of uj at the 
two infinities 




(3.10) 


At the same time from (3.5) we get 


afj — ±Iij + iOij 


lij 



Qi(v)Qj(v) - Qi(v)Qj(v) , 


(3.11) 


which implies that 

= -ilki cot ttMi. (3.12) 

We see that the zero modes can be also computed from the values of Q and Q on the cut. 

Note also that the r.h.s. is not explicitly antisymmetric. Imposing the antisymmetry 
gives 

IkiicotnMi — cotvrMfc) = 0, (3.13) 

so either Iki = 0 or cotvrM; = cotnAdk- As Pf a; = 1, all Iki can not be equal to zero 
simultaneously. Having Iki non-zero implies quantization of charges: for example, the 
choice Ii 2 7 ^ 0 and 134 7 ^ 0 , which is consistent with perturbative data, requires cot vrMi = 
cot 7 rM 2 and cot vrMa = cot and so Si, S 2 have to be integer or half integer. In section 
4 we will see how to relax this condition and do an analytic continuation in the spin S'! to 
the whole complex plane. 


3.3 Step 3: Reducing to an optimization problem 

Having ujij and at hand we can try to impose the remaining equations of the QSC 
(2.17). We notice that there are two different ways of computing Qj, which should give 
the same result when we have a true solution: (2.9) and (2.17). Their difference, which at 
the end should be zero, is 

Fi{u) = P“(u) Qai^iu + i/2) + io.jiu) Q“l^(u + f/2)P,(u) . (3.14) 
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The problem is now to find Ca,n for which Fi{u) is as close as possible to zero. Here we have 
some freedom in how to measure its deviation from zero, but in our implementation we 
use the sum of squares of Fi at the sampling points ua- Then the problem reduces to the 
classical optimization problem of the least squares type. In our implementation we found 
it to be particular efficient to use the Levenberg-Marquardt algorithm (LMA), which we 
briefly describe in the next section. The LMA is known to have a Q-quadratic convergence 
rate, which means that the error Cn decreases with the iteration number n as fast as 
The convergence is indeed so fast that normally it is enough to do 2 or 3 iterations to get 
the result with 10 digits precision. We give some examples in the next section. 

Levenberg-Marquardt algorithm Our problem can be reformulated as follows: given 
a vector function fi{ca) of a set of variables Ca (which we can always assume to be real) 


find the configuration which minimizes 

S(c.) ^ 5] |/.(cj|" . (3.15) 

i 

Assuming we are close to the minimum we can approximate fi by a linear function; 

fi{ca) - fi{Ca) + (Ca “ Ca)Jai , Jai = dc^fi{Ca) (3.16) 

which gives the following approximation for S{ca)'- 

S{Ca) = [fi{Ca) + (Ca “ Ca)Jai] [/i(Ca) + (Ca “ Ca)Jai] (3.17) 

The approximate position of the minimum is then at dc^S = 0 for which we get 

Jai [/i(Ca) T (Ca CQi)TQij] T [fi{Ca) T (Cq Ca)Jai\ Jai — 0 (3.18) 

from which, in matrix notation, 

c = c-{jF + Jj^)-^{Jf + Jf) . (3.19) 


We see that for this method we should also know the derivatives of our Fa{u) w.r.t. the 
parameters Ca^n, which in our implementation we find numerically by shifting a bit the 
corresponding parameter. 

In some cases, when the starting points are far away from the minimum, the above 
procedure may start to diverge. In such cases one can switch to a slower, but more stable, 
gradient descent method for a few steps. The Levenberg-Marquardt algorithm provides a 
nice way to interpolate between the two algorithms by inserting a positive parameter A 
into the above procedure, 

Cn+l =Cn- {JJ^ + JJ^ + AiyHJf + Jf) . (3.20) 

The point is that for large A this is equivalent to the gradient descent method. Thus one 
can try to increase A from its zero value until 5(c,i+i) < S{cn) and only then accept the 
new value Cn+i- This helps a lot to ensure stable convergence. 

In the next section we demonstrate the performance of our numerical method by ap¬ 
plying it to the twist-2 operators in sl(2) sector. 
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3.4 Implementation for the s[(2) Sector and Comparison with Existing Data 

The s[(2) sector in the QSC framework. Although our method can be used for any 
state of the N = A SYM theory, the examples we provide in this paper are for states in the 
s[(2) snbsector. In this section we will discuss the physical operators which have integer 
spin, and demonstrate our numerical method in action for the Konishi operator. Then 
in section 4 we will show how the algorithm works for other states with S no longer an 
integer. 

Let ns sketch the needed information about the QSC in the s[(2) subsector (more 
details can be found in [2, 1, 15]). Operators in this sector have only three non-zero 
quantum numbers: spin S = Si, twist L = Ji and conformal dimension A. Twist-L 
single-trace operators of this subsector can be schematically represented as 

0 = Tt [D^Z^) + ..., (3.21) 


where D is a light-cone derivative, Z is a scalar of the theory and the dots stand for 
permutations. 

For such states the QSC enjoys several simplifications. First, quantities with upper 
and lower indices are now related to each other: indices can be raised or lowered using a 
simple matrix® 


/ 0 

0 

0 

-l\ 

0 

0 

1 

0 

0 

-1 

0 

0 

VI 

0 

0 

0 ^ 


(3.22) 


for example. 




(3.23) 


It is also easy to show that in this sector ojij should satisfy ujn = a ;23 in addition to 
antisymmetry. Second, Q-functions have now definite parity in u, which decreases the 
number of expansion coefficients in series like (3.1) two times. Finally, one can write down 
a simplified version of asymptotics (2.1): 

Pa ~ (Ain-^/2-1^ (3.24) 

_ , A-S A+S-2 A+S -A+S-2 . 

Qa ^ (Biu 2 ,B2U 2 ,B3U 2 ^B^u 2 ). (3.25) 


Coefficients Aa and Bi are related to the global charges L, S, A (see [1, 2]). 

Improved implementation: skipping the computation of cn’s. We have mentioned 
before that the simplification of the QSC achieved in [26] should allow to significantly 
improve the iterative procedure, as one can avoid calculating oJij. Here we present this 
improvement for symmetric states in the si{2) sector. Let us briefly recall the trick used in 
[26] to eliminate w’s. For the states we consider, each of the Pa(^t) functions is either even 
or odd. Then, as follows from the 4th order hnite difference eqation on Qj with coefficients 

^Notice that the choice of x fixes the leading order of Qa\i up to a rescaling. 


- 12 - 








Figure 4: Convergence of the algorithm. The error as measured by the value of (3.15) 
reduces at the quadratic rate e„ ~ e“'^ ^ as a function of the iteration number. In most cases 
our program managed to find the solution from a very remote starting point. On the picture we 
started from all free parameters Ca,n set to zero and with the initial value for the energy Aq = 4.1. 
After 12 iterations it correctly reproduced A = 4.4188599 at A = 167 r^( 0 . 2 )^ ~ 31.6. With each 
iteration taking about 1.5sec the whole procedure took about 20 sec on a Laptop with Intel i7 
2.7GHz processor. 


built from P’s Qi{—u) satisfy the same finite difference equation as Qj(rt). Thus each of 
the former can be expressed as a linear combination of the latter with periodic coefficients: 

Qi{u) = . (3.26) 

We work in the basis where Qj have pure power-like asymptotics at large u, non-coinciding 
for general values of global charges. It is easy to see that at large u in this basis fij(u) 
should be constant and diagonal. At the same time, (3.26) allows us to relate Qi(—u) and 

Qd 

Qi(u) = ajQj(-u), aj = (3-27) 

where y is defined in (3.22). The functions Qi(— m) and Qi(u) have the same analytical 
properties, so a*- should be i-periodic and analytic. One should also take into account that 
only ll>i 2 and oo'n are non-zero at infinity, thus many elements of have to be zero. For 
indices 1 and 3 we finally get the key relations which appear to be sufficient to close the 
QSC equations: 

Qi(tt) = ai3Q3 (-M), Q3(w) = a3iQi(-«) (3.28) 

Consistency of these two equations also implies that 0:13 = 1/031 = a. Note that Qi(m) can 
be constructed in our algorithm as The equation above tells us that it should be 

proportional to Q 3 (—u) with unknown constant of proportionality. This requirement can be 
also phrased as a minimization problem. For that let us evaluate the ratio Qi (^^)/Q3(-«) 

®its explicit form is given e.g. in (3.2) of [15] 
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at sampling points Uk on the cut [—2g, 2g] and compute its variance, 


M 

k=l 

where the constant B is the mean value, 


Qa\i{uk + i/2)P°‘{uk) 
Qs( 


1 ^Qa\i{uk + i/2)P°-{uk) 
^ ^ Q3(-«fc) 


(3.29) 


(3.30) 


On the true solution of the QSC this ratio is a constant so the variance should be zero, i.e. 
5 = 0. Thus our goal is to minimize the function 5, and for this we again use the Levenberg- 
Marquardt algorithm described above. This gives the desired numerical prediction for the 
coefficients Ca,n parameterizing the P-functions. 

The main performance gain stems from the fact that as we do not compute w’s, we no 
longer need to calculate the integrals (3.5) and (3.11). This improved method can also be 
used for non left-right symmetric states as welF, with simple modifications (details of this 
generalization will be presented elsewhere). 


Implementation for Konishi Here we discuss the convergence on a particular example 
of the Konishi operator which corresponds to 5 = 2, L = 2. The reason we start from this 
operator is because it is very well studied both analytically at weak and strong coupling 
and also numerically. So we will have lots of data to compare with. 

To start the iteration process described in the previous sections, we need some rea¬ 
sonably good starting points for the coefficients Ca^n- For the iterative methods, like, for 
instance, Newton’s method, good starting points are normally very important. Depending 
on them the procedure may converge very slowly or even diverge. We made a rather radical 
test of the convergence of our method by setting all coefficients to zero, except the leading 
ones, which are fixed by the charges. For A we took the initial value 4.1 at the value of 
’t Hooft coupling g = 0.2. To our great surprise it took only 12 steps to converge from 
the huge value of S{ca) ~ 10+^ (defined in (3.15)) to 5(ca) ~ 10“®. The whole process 
took about 20 seconds on a usual laptop (see Fig. 4), producing the value A = 4.4188599, 
consistent with the best TEA estimates [18, 24]. 

After that we used the obtained coefficients as starting points for other values of the 
coupling to produce a large volume of data, part of which is shown in Table 1. All the 
values obtained are consistent with the TEA results within the precision of the latter, being 
considerably more accurate at the same time. 

The reason for such an excellent convergence is the Q-quadratic convergence rate of 
the algorithm we use. It means that the number of exact digits doubles with each iteration, 
or that the error decreases as at the step n, if the starting point is close enough. 

What is perhaps surprising is that the algorithm converges from a very remote starting 
point. 

^It has been applied to study the BFKL regime with nonzero conformal spin in [27] 
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vA 

477 

A5=2(A) 

vA 

477 

A5=2(A) 

0.1 

4.115506377945221056840042671851 

0.2 

4.418859880802350962250362876243 

0.3 

4.826948662284842304671283425271 

0.4 

5.271565182595898008221528540034 

0.5 

5.712723424787739030626966875973 

0.6 

6.133862814488691819595425762346 

0.7 

6.531606077852440195886557953690 

0.8 

6.907504206024567515828872789717 

0.9 

7.2641695874391127748396398539 

1 

7.60407071704738848334286555 

1.1 

7.9292942641568451632186264 

1.2 

8.241563441147703542676050 

1.3 

8.54230287229506674486342 

1.4 

8.8326999393163090494514 

1.5 

9.11375404891588560886 

1.6 

9.386314656368554140399 

1.7 

9.65111042653013781471 

1.8 

9.9087717085593508789 

1.9 

10.1598480131615473641 

2 

10.4048217434405061127 

2.1 

10.6441190951617575972 

2.2 

10.878118797537726796 

2.3 

11.107159189584305149 

2.4 

11.331544000504529107 

2.5 

11.551547111042160297 

2.6 

11.76741650605722239 

2.7 

11.97937757952067741 

2.8 

12.18763591669137588 

2.9 

12.3923796509149519 

3 

12.5937814717988565 

3.1 

12.7920003457144898 

3.2 

12.9871829973986392 

3.3 

13.1794651919629055 

3.4 

13.368972849208144 

3.5 

13.555823016292914 

3.6 

13.740124720157966 

3.7 

13.921979717391474 

3.8 

14.101483156227149 

3.9 

14.278724162943763 

4 

14.45378636296056 

4.1 

14.62674834530641 

4.2 

14.79768407780976 

4.3 

14.96666327925592 

4.4 

15.13375175384302 

4.5 

15.29901169250472 

4.6 

15.4625019450274 

4.7 

15.6242782663505 

4.8 

15.7843935399844 

4.9 

15.942897981092 

5 

16.099839321454 


Table 1: Conformal dimension of Konishi operator 


Another indicator of the convergence is the plot of Q computed in two different ways, 
i.e. (2.9) and (2.17). On the true solution of the QSC both should coincide. On Fig. 5 
we show how fast the difference between them vanishes with iterations, i.e. how fast we 
approach the exact solution of the QSC. 

Our numerical data allows to verify with high precision the strong coupling analytic 
prediction for the Konishi dimension obtained in [28, 18, 21, 19, 20, 23, 25, 12], 


A an 

^konishi 




= 2 A^/'^ + 


Ai/4 

2 

AV4 


A3/4 A5/4 

3.106170709 15.48929958 

+ -- + ••• 


(3.31) 


A3/4 


A5/4 


Fitting our data we find 


= -0.0000000475AV4 + + ^ + ^ + ... 


^konishi 


Al/4 


A3/4 A5/4 


(3.32) 
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Figure 5: Q-functions at the first several iterations. Here we show how Q 3 converges to 
the solution in just four iterations when calculating the Konishi anomalous dimension. At each 
picture solid and dashed blue lines show Q 3 slightly below the cut calculated with (2.17) and (2.9) 
respectively, which should coincide on the solution. Red lines show the same slightly above the cut. 

where we see that the first coefficient is reproduced with high precision. Fixing its value 
to match the analytic result, we find that the precision of the fit at subsequent orders in 
1/A improves. Proceeding in this way, we have confirmed all coefficients in the analytic 
prediction (3.31) with at least 3-4 digits precision, which is far more than the TEA data 
allows to get at these high orders in the expansion. 

In the next section we discuss the analytic continuation in S away from its integer 
values. This is an important calculation which bring us to a highly accurate numerical 
estimate for the pomeron intercept at finite coupling — a quantity which can be studied 
exclusively by our methods. 

4. Extension to Non-Integer Lorentz Spin 

In this section we explain which modifications are needed in order to extend our method 
to non-integer values of spin S, and give two specific examples of calculations for such S. 

4.1 Modification of the Algorithm for Non-Integer Spin 

First we need to discuss how the procedure of fixing zero modes of w’s described in section 

3.2 is modified for non-integer S. The main difference stems from the fact that analytic 
continuation to non-integer S changes the asymptotic behavior of coij at large u, as de¬ 
scribed in [ 12 , 15, 30]. While for integer S asymptotics of uj are constant, for non-integer S 
some components of u have to grow exponentially. Without this modification the system 
has no solution: indeed, in section 3.3 we assumed constant asymptotics of all w’s and 
derived quantization condition for global charges. 

A minimal modification would be to allow exponential asymptotics in one of the com¬ 
ponents of UJ. In order to understand which of the components can it be, let us recall the 
Pfaffian constraint satisfied by ujij 

Pf w = wi2u;34 - a;i3ti;24+ 1^14 = 1- (4.1) 

First, it is clear a;i 4 alone can not have exponential asymptotics. Second, in the case 
of integer S both uji 2 of 0:34 are non-zero constants at infinity [15, 12 ]; then shifting S 
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infinitesimally away from an integer we see that it would be impossible to satisfy the 
condition (4.1) if we allow one of them to have exponential asymptotics at infinity: this 
exponent will multiply the constant in the other one. So the only two possibilities left are 
a;i 3 and UJ 24 , which are both zeros at inhnity for integer S. From perturbative data we know 
that it is UJ 24 which should have exponential asymptotics. Thus we formulate the “minimal” 
prescription for analytic continuation of Q-system to non-integer S: asymptotic has 

to be allowed in uj 24 - This prescription was tested thoroughly on a variety of examples 
[31, 32, 15, 12, 30], but it would be interesting to derive it rigorously and generalize it to 
states outside of the s[( 2 ) sector. Of course, one can also consider adding exponents to 
more than one component of cjij: in this case the solution will not be unique. A complete 
classihcation of solutions of Q-system according to exponents in their asymptotics might 
be interesting. For example it is known that allowing for an exponent in some other 
components corresponds to the generalized cusp anomalous dimension [12, 33]. 

Because of the exponential asymptotics of LU 24 , the argument in section 3.2, which fixes 
the zero modes of uj, has to be modihed. First, formula (3.12) still holds true for i = 1 or 
i = 3, as UJ 24 does not enter anywhere in the derivation. Thus 


CO 12 = -ih 2 cot 


7r(5 + A) 
2 


ijJ^4 — — 11^4 cot 


7r(5-A) 
2 


(4.2) 


Consequently, one can use (3.12) for both cuis and wsi, and reproduce the quantization 
condition (3.13) for global charges, which in this case implies that either A = 0 or W 13 = 0. 
Equation (3.12) can also be used for t<;i 4 and W 23 (which are equal) and imposes that either 
A = 0 or wi 4 = 0. 

It remains to hx the zero mode in u; 24 , for which we use an ansatz 

a;24 = + 02 + 036 “^™. (4.3) 


The constants 04 , 02,03 can be found from the Pfaffian constraint (4.1). To this end we 
expand the hyperbolic cotangent in (3.5) to get 

Wij = + lij + 2e + 2e + ..., u —>■ + 00 , (4-4) 

where the terms of the expansion are integrals similar to lij with additional factors of 
^ 2 ttu qj, g 47 ru integrand^. Analogous expansion can be obtained at n —)• — 00 . Then 

plugging these expansions into (4.1) we get formulas for the coefficients 04 , 02 , 03 . For 
example, 

l + Zi^ ("1 +jcot (l-zcot^^^^) 

04 = 2 i -^-T . ( 4 . 5 ) 

''13 

With these modifications we can reconstruct all ujij including the zero modes and then 
proceed with our algorithm as in the case of integer S. 


4.2 Exploring Complex Spin 

In this section we briefly describe the results of our numerical exploration of A (S') as an 
analytic function of a complexified spin S. As explained in the previous section the general¬ 
ization of our numerical method to arbitrary values of spin requires minimal modifications 
® Actually, these integrals can be evaluated analytically in terms of Bessel functions 
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Figure 6 : Section of the Riemann surface S{A) along Im A = 0 for different values of coupling g. 
The upper two solid curves, shown in black and grey, represent the well-known BFKL eigenvalue 
as a function of A, whereas the lower two come from the unphysical sheet which can be accessed 
from the upper one by going through the cuts. The dashed line shows the zero-coupling limit of 
the curve. Orange dots mark BPS states Tr(ZZ). 


of our main code. Thus we are able to generate numerous values of the anomalous di¬ 
mension for any S with high precision in seconds. In fact both S and A enter the QSC 
formalism on almost equal footing and we can also switch quite easily to finding S for given 
A. This is what is adopted in the vast literature on the subject and what we are going 
to consider below. This viewpoint is particular convenient due to the symmetry A —>■ — A 
which makes the pictures particularly nice. 

Starting from S = 2 (Konishi operator) we decreased the value of S' or A in small 
steps using the solution at the previous step as a starting point for the next value. In 
this way we built the upper two curves on Fig. 6 . Let us point out one curious technical 
problem - one can see for instance from (4.5) that the lines S = ±A -|- Z are potentially 
dangerous due to the divergence. In fact we found that near these dangerous points on the 
line the factor / 12 L 34 also vanishes canceling the potential divergence. This however affect 
the convergence “radius” of our iterative procedure and we found it quite complicated to 
cross those lines, even though in very small steps we were able to reach close to them. The 
way out is to go around these lines in the complex plane A. To make sure there is no true 
singularity or branch point we also explored a big patch of the complex plane A, indeed 
finding some branch points, but deep inside the complex plane, having nothing to do with 
these lines. For example when g = 0.2 we found 4 closest branch points at roughly ±1 ± i, 
see Fig. 1. By making an analytic continuation (described above) through those cuts we 
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Figure 7: The BFKL intercept j as a function of coupling A. The red solid line with tiny red dots 
is obtained by our numerical procedure. It interpolates perfectly between the known perturbative 
predictions (the blue dashed lines) at weak [34, 35] and strong coupling [36, 37, 38, 12]. 


found another sheet of the Riemann surface ^(A). On this sheet we have found four cuts: 
two are connecting it to the first sheet and two other ones, located symmetrically on the 
imaginary axis, lead to further sheets. We expect an infinite set of sheets hidden below 
and also more cuts on both sheets outside of the area that we have explored. 

It is instructive to see how this Riemann surface behaves as g ^ 0. First, the real 
parts of branch points on the physical sheet are very close to ±1, but the imaginary part 
goes to zero. Thus at inhnitely small g the cuts collide, isolating the region [Re A| < 1 
from the rest of the complex plane. These two separated regions become then the areas 
of applicability of two different approximations: for [Re A| > 1 one can apply the usual 
perturbation theory and Beisert-Eden-Staudacher Asymptotic Bethe Ansatz, whereas the 
region [Re A| < 1 is described by BFKL approximation and so-called Asymptotic BFKL 
Ansatz [30]. 

The presence of the cut can be to some extent deduced from perturbative perspective 
in each region: in the regime of usual perturbation theory 

A{S) = 2 + S-8g^Hs + Oig^), (4.6) 


where Hg is the harmonic number. It has poles for all negative integer values of 5 — 
these poles are weak-coupling remnants of the cuts we see at finite coupling. In the BFKL 
regime one should instead look at the leading order BFKL equation [39, 40, 34] 


5(A) = -1 + 45^ 


V’ 


1 + A 


-b ijj 


I - A 


-2V^(1) 




(4.7) 
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To make sense of this equation one has to take the limit g ^ 0 , S' —>• — 1 so that the l.h.s 
stays finite. Then the '^/’-functions in the r.h.s generate poles at odd values of A, which, 
again, are cuts degenerated at weak coupling. 

Fig 6 represents a section of the Riemann surface by the plane Im u = 0, i.e. de¬ 
pendence of S on A for real A, which, of course, consists of two curves, originating from 
the two sheets we explored. At weak coupling the upper curve becomes piecewise linear, 
approaching different parts of the dotted line: for |A| > 1 it coincides with S = ±A — 2 
and for |A| < 1 it becomes S = — 1. One could expect a similar piecewise linear behavior 
for the lower curve: it approaches S = ±A — 2 for |A| < 1, approaches S' = 0 in some 
region outside of |A| < 1 and becomes a certain linear function even further away from 
A = 0. It would be interesting to explore the complete analytic structure of this Riemann 
surface, and understand what describes its asymptotics when g ^ 0. It should produce a 
hierarchy of “Asymptotic Bethe Ansatze” each responsible for its own linear part of the 
limiting surface. 

4.3 BFKL Pomeron Intercept 

The pomeron intercept j(A) is a quantity which relates spectrum of single-trace operators 
and scaling of high energy scattering amplitudes in the Regge regime [41]. This regime is 
particularly interesting, since it establishes a connection between results in AA = 4 SYM 
and multicolor QCD: the non-trivial leading order of so-called BFKL eigenvalue is the same 
in two theories, and in the higher orders A/" = 4 SYM is expected to reproduce at least the 
most complicated part of the QCD result. 

Our goal is to demonstrate the universal power of our approach by giving an extremely 
precise numerical estimate for this important quantity at finite coupling in a wide range of 
couplings. 

One defines the intercept as j = ^(A = 0) -|- 2, where S is the spin of the twist-2 
operator such that A (S') = 0. Having formulated the problem like this, we can in principle 
apply the algorithm described in section 3 to find the correct value of S, while keeping 
A at zero. However, one may already suspect that the point A = 0 is special. Indeed, 
we know that for any solution of QSC there is always another one related by A —)• —A 
symmetry. At the level of Qj functions this allows simultaneously interchanging Qi -H- Q3 
and Q20Q4 as one can see from the asymptotics ( 3 . 25 ). From this we see that at small A 
two different solutions of QSC (related by the symmetry) approach each other, making the 
convergence slow, exactly like Newton’s method becomes inefficient for degenerate zeros. 
In other words, in the limit A —>• 0 the Q’s related by the symmetry become linearly 
dependent in the leading order. Furthermore, since the matrix Qa\i should stay invertible, 
the leading coefficients Bi of asymptotic expansion of Qj diverge at A —>• 0. 

The way out is to perform a linear transformation of Q’s preserving the equations: 
it will replace two of them by linear combinations Q 3 — 7Q1 and Q 4 -|- 7Q2 with some 
coefficient 7, so that the divergent leading order cancels and the four functions Qi become 
linearly independent. 

For the gauge choice in which Bi = B2 = 1 the transformation acts on i-indices of 
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Q-functions with a matrix®. 


/ 1 0 0 o\ 

0 100 _ _ z(5-4)(5-2)5(5 + 2) 

-7 0 1 0 ’ ^ “ 16(5- 1)2A 

\ 0 7 0 1 / 


(4.8) 


One can check that rotation by this matrix will render Qa\i finite and linearly inde¬ 
pendent, and moreover, preserve relations (3.23). After this one can apply the standard 
procedure from section 3 with the only modification that the large u expansion of Qa\i will 
contain logu/rt"' terms in addition to the usual Xju^. 

Having done this, we can readily generate lots of numerical results. In particular we 
built numerically the function j(A) which interpolates perfectly between the weak and 
strong coupling predictions. We have found j{X) with high precision (up to 20 digits) for 
a wide range of’t Hooft coupling (going up to A ~ 1000). The results are also summarized 
in the Table 2. This table represents a small portion of all data we generated, which is 
available by request. 


a/A 

477 

J(A) 

vA 

477 

J(A) 

0. 

1.000 000 000 000 000 000 0 

0.1 

1.101144 978 997 772 874 8 

0.2 

1.301794032 258 782 208 7 

0.3 

1.470 445 240 989187630 6 

0.4 

1.587128 066 254129 730 4 

0.5 

1.666 438 709 974061852 3 

0.6 

1.721917 842 815 631353 9 

0.7 

1.762 239 296 816 453 814 3 

0.8 

1.792 626 253 069403 59 

0.9 

1.816 252 952 807 28411 

1. 

1.835109 464032173 0 

1.1 

1.850 489 553 739 522 8 

1.2 

1.863 264 346 392 6404 

1.3 

1.874039 320 799 460 

1.4 

1.883 247 290 966 33 

1.5 

1.891205 346 040 23 

1.6 

1.898150 851852 49 

1.7 

1.904 264 892 92817 

1.8 

1.909 687 948 27174 

1.9 

1.914 530 628 017 38 

2. 

1.918881 1873049 

2.1 

1.922 810 887 750 

2.2 

1.926 377 890 67 

2.3 

1.929 630129 41 

2.4 

1.932 6074591 

2.5 

1.935 343 2872 


Table 2: Numerical data for the pomeron intercept for various values of the ’t Hooft coupling. 
All digits are expected to be significant. 


At strong coupling we can confirm the analytic predictions obtained in [36, 37, 38, 12], 


111 6C(3)+2 18C(3) + ^f 39C(3) + ^ 

A 4 A3/2 A2 ^ A'5/2 ^ A3 ^ 

1 ^ 0.2500000000 ^ 9.212341419 ^ 27.27764926 , 62.84896922 

A+ + A2 ^ A572 ^ A3 ^ 


®This is a particular case of H-transformations described in section 4.1.3 of [2] 


= 2 - 


= 2 - 


2 

AV2 

2 

AV2 
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Fitting our numerical results, we find 


juum _ jan _ 


-0.000000687 + 


0.000149 0.0146 0.854 


AV2 


+ 


A 


+ 


A3/2 


33.105 892.72 

A2 ^ A5/2 


(4.10) 

17093.5 
A3 ■ 


As for the Konishi anomalous dimension (see (3.32)), we see that the leading coefficient 
is reproduced with high precision, while at the next orders the precision decreases and 
we get coefficients very different from the analytic prediction. However, fixing the leading 
coefficient to match the analytic result, we again found that the precision of our fit at 
the next orders in 1/A increases. Gradually reproducing the coefficients in this way, we 
confirmed all coefficients in the analytic result (4.9) with the precision of at least 3-4 digits, 
except for the A3 term for which we could not get a stable fit with the data we have so far. 

Moreover, we generated ~ 100 points with small g in the range 0.1... 0.017 each with 
more than 20 digits precision. Fitting this data with powers of we found 


j = 1 + 11.09035488895912/ - 84.0785668075/ - 2543.0481652/ 156244.8086c/® 


where the first 3 terms are known analytically from Feynman diagram perturbation theory 
calculations [34, 35] and their numerical values coincide in all digits with our prediction 
above. The last two terms give our numerical prediction for the numerical values of the 
NNLO and NNNLO BFKL pomeron intercept. Our fit also gives predictions for the higher 
corrections but with a smaller precision. 


5. Conclusions and Future Directions 

In this paper we have demonstrated that in addition to their analytic power, the QSC 
equations can give highly precise numerical results at finite coupling. We develop a numer¬ 
ical procedure which applies to generic single trace operators and as such it is unique in 
its kind. Furthermore, the algorithm converges at a remarkably high rate which gives us 
access to high numerical precision results - up to 20 digits or even more in a few iterations. 

The efficiency of our method is demonstrated on the example of s[(2) sector opera¬ 
tors. We also formulated how to extend our procedure to non-integer quantum numbers. 
We studied the twist-2 operators for complex values of the spin discovering a fascinating 
Riemann surface (see Fig.l). In addition we reformulated our equations to be directly 
applicable to the BFKL pomeron intercept and evaluated the intercept j with high preci¬ 
sion of up to 20 significant figures. By fitting our data we also gave a prediction for the 
perturbation theory expansion 

j(A) = 1 + 0.07023049277268284 A - 0.00337167607361 A^ (5.1) 

- 0.00064579607573 A® -k 0.0002512619258 A"^ -k ... 

reproducing correctly the first two nontrivial orders [34, 35] and giving a prediction for 
higher orders. 

The range of possible applications of our method is vast. First, it is not limited solely 
to the s[(2) sector of AA = 4 SYM, but is directly applicable to any single trace operators 
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of the theory. It would be interesting to do an explicit example of a numerical calculation 
with our algorithm outside of the s[(2) sector. For example, the wider class of s[(2,C) 
operators (identified in [30]), also exhibiting a BFKL regime, could be a good candidate to 
begin with. Second, we expect our method to be applicable for such non-local operators as 
the generalized cusp anomalous dimension and quark-anti-quark potential, DD-brane and 
other boundary problems [42, 43, 11, 44, 33]. Third, it may be interesting to generalize 
our method to ABJM theory as well as to various integrable deformations of AA = 4 SYM 
theory. 

The numerical results could also be helpful for the analytical exploration of the spec¬ 
trum - for instance, in such regimes as BFKL and at strong coupling, which remains almost 
unexplored, and various limiting cases of the generalized cusp. Furthermore, studying nu¬ 
merical results and the behavior of the generated Q-functions in various limits can reveal 
new analytically solvable regimes. 

To facilitate further applications and development of our method we attach to this 
paper a user-friendly version of our code as a Mathematica notebook. It provides an 
implementation of our algorithm in the simplest case. It may be also useful to convert our 
Mathematica code into a faster language such as C-|—|- or a similar lower-level language 
It should not be difficult as our algorithm is quite simple and only uses some basic matrix 
operations. 
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